library(data.table)
library(igraph)
library(tidyverse)
library(stargazer)
library(tibble)
library(RColorBrewer)
library(digest)
library(viridis)
library(Cairo)  # for saving unicode chars in PDF

set.seed(1)

plot_theme <- function(base_size = 16) {
  theme_set(theme_bw(base_size = base_size))
  theme_update(panel.grid = element_blank(),
               panel.border = element_rect(size = .5),
               panel.grid.major = element_line(size = .1, color = "#f0f0f0"))
}
plot_theme(16)

current_prognosis_label <- "Current Prognosis (One-Year Mortality Risk)"
initial_prognosis_label <- "Initial Prognosis (One-Year Mortality Risk)"

event_type_colors <-
  rev(
    c(brewer.pal(6, "Greys")[1],
      brewer.pal(6, "Blues")[2],
      brewer.pal(6, "Purples")[3],
      brewer.pal(6, "Greens")[4],
      brewer.pal(6, "Oranges")[5],
      brewer.pal(6, "Reds")[6])
  )

hist_separator <- "->"

# import actual history phat and cost --------------------------------------------


hist_data_raw <- fread("event_tree.csv") 


rename_events <- function(x) {
  recode_factor(x, 
                  `High` = "Admission-High",
                  `Low` = "Admission-Low",
                  `drugs` = "Drug Therapy",
                  `radiotherapy` = "Radiation Therapy",
                  `ER` = "ED Visit",
                  `start` = "Initial Diagnosis"
                  )
}

hist_data <- hist_data_raw %>% 
  rename(hist_id = V1, 
         cancer_type = CNRS_topo_main_groups,
         mean_phat0 = phat0_avg,
         mean_phat_current = phat_avg,
         mean_avg_monthly_fwd_cost = fwd_cost_avg,
         mean_avg_monthly_fwd_cost_excl_index = fwd_cost_no_current_avg,
         history_depth = EVE_event_number)  %>% 
  mutate_at(.vars = vars(starts_with("EVE_")), .funs = rename_events) %>% 
  mutate(hist_id = paste0("H", hist_id),
         hist_name = paste(cancer_type,
           EVE_lag7_category_new,
           EVE_lag6_category_new,
           EVE_lag5_category_new,
           EVE_lag4_category_new,
           EVE_lag3_category_new,
           EVE_lag2_category_new,
           EVE_lag1_category_new,
           EVE_category,
           sep = hist_separator
         ) %>% 
           gsub(pattern = paste0("NA", hist_separator), replacement = "", fixed = T, x = .),
         current_event_type = rename_events(EVE_category))
all_cancer_types <- unique(hist_data$cancer_type)
main_cancer_types <- c(Breast = "Breast",
                       `Prostate gland` = "Prostate",
                       Colon = "Colon",
                       Brain = "Brain",
                       `Bronchus and lung` = "Lung",
                       Pancreas = "Pancreas")
actions <- unique(hist_data$current_event_type) 
                

#mold histories into a tree ----------------------------------------------


parent_to_child <- function(parent_histories, actions, depth = 1) {
  ans <- data.frame()
  while (depth > 0) {
    parents <- rep(parent_histories, each = length(actions))
    children <- paste(parents, actions, sep = hist_separator)
    ans <- bind_rows(ans,
                     data.frame(parent = parents,
                                child = children))
    parent_histories <- children
    depth <- depth - 1L
  }
  return(ans)
}

full_tree_edges <- parent_to_child(all_cancer_types,
                                   actions, depth = 7)


# bind edges to get a crosswalk from parents to children histories
hist_list <- hist_data %>% select(hist_id, hist_name)
hist_edges <-
  rename(hist_list, hist_id_parent = hist_id, parent = hist_name) %>%
  inner_join(full_tree_edges) %>%
  inner_join(rename(hist_list, hist_id_child = hist_id, child = hist_name)) %>%
  select(hist_id_parent, hist_id_child)

# note orphan stubs
orphan_hist <- setdiff(unique(hist_data$hist_id), unique(with(hist_edges, union(hist_id_parent, hist_id_child))))
hist_data %>% filter(hist_id %in% orphan_hist) %>% pull(hist_name)


# merge histories with edges
hist_data_parent_child <- rename(hist_data, hist_id_child = hist_id) %>%
  left_join(hist_edges, by = "hist_id_child") %>%
  left_join(rename(hist_data, hist_id_parent = hist_id), by = "hist_id_parent",
            suffix = c("_child", "_parent"))

fwrite(hist_data_parent_child, "hist_parent_child.csv")
hist_data_parent_child <- fread("hist_parent_child.csv")


# history descriptives ----------------------------------------------------

hist_desc <- hist_data %>% group_by(history_depth, current_event_type) %>% summarise(histories = n(), patients = sum(N)) 

#fig 2a
hist_desc %>% 
  ggplot(aes(x=as.factor(history_depth), y = patients, group = current_event_type, fill = current_event_type)) +
  geom_col(color = "gray80", size = .25) + 
  labs(x = "Event Number", y = "Number of Cases", fill = "Event Type") +
  theme(text = element_text(size = 21)) +
  scale_fill_manual(values = event_type_colors) 
ggsave("ev_fig_number_of_patient_by_event_type_within_history_depth_stacked.pdf", w= 12, h =6 )

#fig 2b
hist_desc %>% ggplot(aes(x=as.factor(history_depth), y = patients, group = current_event_type, fill = current_event_type)) + 
  scale_fill_manual(values = event_type_colors) +
  theme(text = element_text(size = 21)) +
  geom_bar(stat = "identity", position = "fill", color = "gray80") + 
  labs(x = "Event Number", y = "Share of Cases",  fill = "Event Type") +
ggsave("ev_fig_number_of_patient_by_event_type_within_history_depth_percent.pdf", w= 12, h = 6)


hist_data[cancer_type %in% names(main_cancer_types), sum(N)]

hist_data %>% 
  filter(cancer_type %in% names(main_cancer_types)) %>% 
  mutate(cancer_type = recode_factor(cancer_type, !!!main_cancer_types)) %>%
  group_by(history_depth, current_event_type, cancer_type) %>% 
  summarise(histories = n(), patients = sum(N)) %>% 
  ggplot(aes(x=history_depth, y = patients, group = current_event_type, fill = current_event_type)) + 
  geom_bar(stat = "identity", position = "fill") + 
  labs(x = "Event Number", y = "Share of Cases", fill = "Event Type") +
  facet_wrap(cancer_type ~.) +
  scale_fill_manual(values = event_type_colors) +
  ggsave("ev_fig_number_of_patient_by_event_type_within_history_depth_and_main_cancer_type_percent.pdf", w= 10, h = 7.5)


# risk/spending dynamics -----------------------------------------------------------

# binned version
binned_current_previous <- hist_data_parent_child %>%
  filter(history_depth_child >= 1) %>%
  group_by(current_event_type_child) %>% 
  mutate(mean_phat_current_parent_bin = ntile(mean_phat_current_parent, 10)) %>% 
  group_by(mean_phat_current_parent_bin, current_event_type_child) %>% 
  summarize(mean_phat_current_child_binned = weighted.mean(
       x = mean_phat_current_child),
       mean_phat_current_parent_binned = weighted.mean(
         x = mean_phat_current_parent),
       N_child = sum(N_child)) 


binned_current_previous %>%  ggplot(aes(
    x = mean_phat_current_parent_binned,
    y = mean_phat_current_child_binned,
    color = current_event_type_child,
    shape = current_event_type_child
  )) +
  geom_abline(slope = 1, intercept = 0, lty = "dashed") +
  geom_point(size = 5) + 
  geom_smooth(method = "lm", se = F) +
  labs(x = "Prognosis (One-Year Mortality Risk)\nat the Start of Previous Event",
       y = "Prognosis (One-Year Mortality Risk)\nat the Start of Current Event",
       shape = "Current\nEvent",
       color = "Current\nEvent") +
  scale_x_continuous(limits = c(0,1), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  scale_color_manual(values = event_type_colors) +
  theme(panel.grid.minor = element_blank())
ggsave("ev_fig_phat_parent_phat_child_by_event_type_binned.pdf", w = 10, h = 7.5)

# sample size counts:
binned_current_previous %>% ungroup %>% 
  group_by(current_event = current_event_type_child) %>% 
  summarise(n = sum(N_child)) %>% fwrite(file = "ev_fig_phat_parent_phat_child_by_event_type_binned_sample_sizes.csv")
  
  

# deltas ------------------------------------------------------------------

sum_diffs <- hist_data_parent_child %>%  
  mutate(phat_diff = mean_phat_current_child - mean_phat_current_parent,
         fwd_cost_excl_diff = mean_avg_monthly_fwd_cost_excl_index_child - 
           mean_avg_monthly_fwd_cost_excl_index_parent) %>%
  select(cancer_type = cancer_type_parent, 
         phat_diff,
         fwd_cost_excl_diff,
         current_phat = mean_phat_current_child,
         N_parent) %>% 
  na.omit %>% 
  group_by(cancer_type) %>% 
  mutate(p5_mean_phat_current_child = cut(current_phat, breaks = c(-Inf,seq(0,1,.2),Inf)),
         q5_mean_phat_current_child = factor(ntile(current_phat, 5), 
                                             levels = 1:5, labels = c("1 (low)", "2", "3", "4", "5 (high)"))) %>% 
  ungroup 
  
baseline_diff_plot <- sum_diffs  %>%  ggplot(aes(x = phat_diff, y = fwd_cost_excl_diff)) +
  labs(x = "Change in Current One-Year Mortality Prognosis",
       y = "Change in Average Monthly Spending,\nOne Year Forward (NIS)") +
  coord_cartesian(ylim = c(-1e4,1e4)) 


baseline_diff_plot + 
  stat_summary_bin(size = .5, alpha=  .7, color = "black", fun = mean, fun.data=mean_se,
                   breaks = quantile(sum_diffs$phat_diff, probs = seq(0, 1, by = 0.05))) + 
  geom_hline(yintercept = 0, color = "gray60", size = .5) +
  geom_vline(xintercept = 0, color = "gray60", size = .5) +
  coord_cartesian(ylim = c(-2e3,2e3), xlim = c(-0.25, .25)) + 
  theme(panel.grid.minor  = element_blank(),
        axis.title = element_text(size = 18)) +
  annotate("text", x = -0.2, y = -2000, color = "gray20", label = sprintf('\u2190 Improving Prognosis'), size = 5) +
  annotate("text", x = +0.2, y = -2000, color = "gray20", label = sprintf('Worsening Prognosis \u2192'), size = 5) 
ggsave("ev_fig_delta_monthly_avg_1yr_fwd_cost_excl_index_over_delta_phat_binned_errorbars.pdf", w = 10, h = 7.5, device = cairo_pdf); dev.off()



ann_text <- data.frame(phat_diff = c(-0.15,+0.15),fwd_cost_excl_diff = -3000,
                       p5_mean_phat_current_child = factor("(0,0.2]\nBest"), 
                       label = c(sprintf('Improving\n\u2190'), sprintf('Worsening\n\u2192')))
sum_diffs %>% 
   mutate(p5_mean_phat_current_child = plyr::revalue(p5_mean_phat_current_child, c("(0,0.2]"="(0,0.2]\nBest", "(0.2,0.4]"="(0.2,0.4]\n", "(0.4,0.6]"="(0.4,0.6]\n", "(0.6,0.8]" = "(0.6,0.8]\n", "(0.8,1]"="(0.8,1]\nWorst"))) %>% 
  group_by(p5_mean_phat_current_child,
           p5_mean_phat_current_child_quantiles = ntile(phat_diff, 10)) %>% 
  summarise(phat_diff = mean(phat_diff), fwd_cost_excl_diff = mean(fwd_cost_excl_diff)) %>% 
  ggplot(aes(x = phat_diff, y = fwd_cost_excl_diff)) +
  labs(x = "Change in Current One-Year Mortality Prognosis",
       y = "Change in Average Monthly Spending,\nOne Year Forward (NIS)",
       subtitle = "Range of current prognosis") +
  geom_hline(yintercept = 0, color = "gray60", size = .5) +
  geom_vline(xintercept = 0, color = "gray60", size = .5) +
  geom_point(size = 2, color = "gray20") +
  geom_smooth(se = F, method = "lm",
              size = 1, color = "black")  +
  coord_cartesian(ylim = c(-3e3,3e3), xlim = c(-0.25, .25)) + 
  theme(panel.grid.minor  = element_blank()) +
  facet_grid(~p5_mean_phat_current_child) + 
  theme(legend.position = "none", axis.text = element_text(size = rel(.6)), panel.grid = element_blank()) +
  scale_x_continuous(breaks = c(-.2,0,.2)) +
  theme(panel.grid.minor  = element_blank(),
        axis.title = element_text(size = 18),
        strip.text = element_text(size = rel(1))) + 
  geom_text(data = ann_text, mapping = aes(label = label), size = 3.5, color = "gray20")
ggsave("ev_fig_delta_monthly_avg_1yr_fwd_cost_excl_index_over_delta_phat_by_current_phat_binned.pdf",
       w = 10, h = 6, device = cairo_pdf); dev.off()



sum_diffs %>% 
  group_by(q5_mean_phat_current_child,
           q5_mean_phat_current_child_quantiles = ntile(phat_diff, 10)) %>% 
  summarise(phat_diff = mean(phat_diff), fwd_cost_excl_diff = mean(fwd_cost_excl_diff)) %>% 
  ggplot(aes(x = phat_diff, y = fwd_cost_excl_diff)) +
  labs(x = "Change in Current One-Year Mortality Prognosis",
       y = "Change in Average Monthly Spending,\nOne Year Forward (NIS)",
       subtitle = "Within-cancer-type quintile of current prognosis") +
  geom_hline(yintercept = 0, color = "gray60", size = .5) +
  geom_vline(xintercept = 0, color = "gray60", size = .5) +
  geom_point(size = 2, color = "gray20") +
  geom_smooth(se = F, method = "lm",
              size = 1, color = "black")  +
  coord_cartesian(ylim = c(-3e3,3e3), xlim = c(-0.25, .25)) + 
  theme(panel.grid.minor  = element_blank()) +
  facet_grid(~q5_mean_phat_current_child) + 
  theme(legend.position = "none", axis.text = element_text(size = rel(.6)), panel.grid = element_blank()) +
  scale_x_continuous(breaks = c(-.2,0,.2)) +
  theme(panel.grid.minor  = element_blank(),
        axis.ticks = element_blank())
ggsave("ev_fig_delta_monthly_avg_1yr_fwd_cost_excl_index_over_delta_phat_by_current_phat_quintile_binned.pdf", w = 10, h = 6)



# age-related exhibhits ---------------------------------------------------


age_group_colors <- c("#4dac26", "#b8e186") 

# import actual history phat and cost --------------------------------------------
rename_events <- function(x) {
  recode_factor(x, 
                `High` = "Admission-High",
                `Low` = "Admission-Low",
                `drugs` = "Drug Therapy",
                `radiotherapy` = "Radiation",
                `ER` = "ED Visit",
                `start` = "Initial Diagnosis"
  )
}


hist_data_raw_age <- fread("event_tree_median_age.csv") 


hist_data_raw_age <- hist_data_raw_age %>% 
  mutate(age = recode_factor(above_median_age,
                             `Below` = "Young",
                             `Above` = "Old",
                             `All` = "All"),
         fwd_cost_current = fwd_cost_avg - fwd_cost_no_current_avg) %>% 
  rename(cancer_type = CNRS_topo_main_groups)



hist_data_age_wide <- hist_data_raw_age %>% 
  pivot_wider(id_cols = c("cancer_type", starts_with("EVE_")),
              names_from = "age",
              values_from = c("N",
                              "phat_avg",
                              "fwd_cost_no_current_avg",
                              "fwd_cost_current")) 

event_by_age_summary <- hist_data_raw_age %>%   group_by(event_number = EVE_event_number, age) %>% 
  summarize(cases = sum(N), current_risk = weighted.mean(phat_avg, w = N), fwd_spending = mean(fwd_cost_no_current_avg, w = N))


event_by_age_summary %>% pivot_wider(id_cols = "event_number", names_from = "age", 
                                     values_from = c("cases", "current_risk", "fwd_spending")) %>% 
  write_csv("ev_AGE_tab_summary.csv")



hist_data_raw_age %>% filter(EVE_category != "start") %>% 
  mutate(event_type = rename_events(EVE_category)) %>% 
  ggplot(aes(x = phat_avg, y = fwd_cost_no_current_avg, color = age, shape = age, size = N)) +
  scale_color_manual(values = age_group_colors)  + 
  stat_summary_bin(aes(color = age), fun='mean', 
                   breaks = quantile(hist_data_raw_age$phat_avg,
                                     probs = seq(0, 1, by = 0.1)), geom='point', size = 2.5) +
  stat_summary_bin(aes(color = age), fun='mean', 
                   breaks = quantile(hist_data_raw_age$phat_avg,
                                     probs = seq(0, 1, by = 0.1)), geom='line', size = .5) +
  labs(x= current_prognosis_label, y = "Average Monthly Spending (NIS)", color = "Patient\nAge", shape = "Patient\nAge") +
  facet_grid(.~event_type, margins = T) + coord_cartesian(ylim=c(0, 2.5e4)) +
  scale_x_continuous(breaks = seq(0, 1, by = .5)) + 
  theme(legend.position = "right",
        text = element_text(size = 18)) +
  ggsave("ev_AGE_fig_fwd_spending_on_current_phat_by_age_and_event_type.pdf", w = 12, h =6)

hist_data_raw_age %>% filter(EVE_category != "start") %>% 
  mutate(event_type = rename_events(EVE_category)) %>% 
  ggplot(aes(x = phat_avg, y = fwd_cost_current, color = age, shape = age, size = N)) +
  scale_color_manual(values = age_group_colors)  + 
  stat_summary_bin(aes(color = age), fun='mean', 
                   breaks = quantile(hist_data_raw_age$phat_avg,
                                     probs = seq(0, 1, by = 0.1)), geom='point', size = 2.5) +
  stat_summary_bin(aes(color = age), fun='mean', 
                   breaks = quantile(hist_data_raw_age$phat_avg,
                                     probs = seq(0, 1, by = 0.1)), geom='line', size = .5) +
  labs(x= current_prognosis_label, y ="Average Monthly Spending (NIS)", color = "Patient\nAge", shape = "Patient\nAge") +
  facet_grid(.~event_type, margins = T) + coord_cartesian(ylim=c(0, 2.5e4)) +
  scale_x_continuous(breaks = seq(0, 1, by = .5)) + theme(panel.grid = element_blank()) +
  theme(legend.position = "right",
        text = element_text(size = 18)) +
ggsave("ev_AGE_fig_index_spending_on_current_phat_by_age_and_event_type.pdf", w = 12, h =6)


hist_data_raw_age %>% rename(history_depth = EVE_event_number) %>% 
  mutate(current_event_type = rename_events(EVE_category)) %>% 
  group_by(age, history_depth, current_event_type) %>% summarise(histories = n(), patients = sum(N)) %>% 
  ggplot(aes(x=history_depth, y = patients, group = current_event_type, fill = current_event_type)) + 
  geom_bar(stat = "identity", position = "fill", size = .5, color = "gray80") +
  scale_fill_manual(values = event_type_colors) +
  facet_grid(. ~ age) +
  labs(x="Event Number", y = "Share of Cases", fill = "Event Type") +
  theme(panel.grid.major = element_blank(),
        text = element_text(size = 16))
ggsave("ev_fig_number_of_patient_by_event_type_within_history_depth_percent_by_age.pdf", w= 10, h = 7.5)

# poor man's reweighting --------------------------------------------------

# calc weights
hist_data_age_compact <- 
  hist_data_raw_age %>% 
  rowwise() %>% 
  mutate(hist_id = digest(algo = "md5", 
                          paste(cancer_type,
                                EVE_lag7_category_new,
                                EVE_lag6_category_new,
                                EVE_lag5_category_new,
                                EVE_lag4_category_new,
                                EVE_lag3_category_new,
                                EVE_lag2_category_new,
                                EVE_lag1_category_new,
                                EVE_category))) %>% 
  select(
    hist_id,
    age = above_median_age,
    cancer_type,
    event_type = EVE_category,
    event_number = EVE_event_number,
    cost_fwd = fwd_cost_no_current_avg,
    cost_current = fwd_cost_current,
    phat_avg,
    N) 

# quick descriptives
hist_data_age_compact %>% group_by(event_type) %>% summarise(mean(phat_avg)) 

# reweighting
tab_reweighted_age <- hist_data_age_compact %>% filter(event_number > 0) %>%  
  group_by(event_type, age, phat_bin = cut(phat_avg, breaks = seq(0,1,.1), include.lowest = T)) %>% 
  summarise(cost_fwd = weighted.mean(cost_fwd, w= N),
            cost_current = weighted.mean(cost_current, w= N),
            N = sum(N)) %>% 
  ungroup %>% 
  pivot_wider(
    id_cols = c("event_type", "phat_bin"),
    values_from = c(starts_with("cost"), "N"),
    names_from = "age"
  ) %>% 
  group_by(`Current Event Type` = event_type) %>% 
  summarise(`Forward Cost, Old` = weighted.mean(cost_fwd_Above, w = N_Above, na.rm = T),
            `Forward Cost, Young` = weighted.mean(cost_fwd_Below, w = N_Above, na.rm = T),
            `Current Event, Old` = weighted.mean(cost_current_Above, w = N_Above, na.rm = T),
            `Current Event, Young` = weighted.mean(cost_current_Below, w = N_Above, na.rm = T),
  )

tab_reweighted_age %>%mutate_if(.predicate = is.numeric, .funs = ~format(., digits = 0, big.mark = ",", scientific = F)) %>% 
  stargazer(summary = F, out = "tab_rewegihted_age.tex", rownames = F)

tab_reweighted_age_total <-
  hist_data_age_compact %>% filter(event_number > 0) %>%  
  group_by(age, phat_bin = cut(phat_avg, breaks = seq(0,1,.1), include.lowest = T)) %>% 
  summarise(cost_fwd = weighted.mean(cost_fwd, w= N),
            cost_current = weighted.mean(cost_current, w= N),
            N = sum(N)) %>% 
  ungroup %>% 
  pivot_wider(
    id_cols = c("phat_bin"),
    values_from = c(starts_with("cost"), "N"),
    names_from = "age"
  ) %>% 
  summarise(`Current Event Type` = "All",
            `Forward Cost, Old` = weighted.mean(cost_fwd_Above, w = N_Above, na.rm = T),
            `Forward Cost, Young` = weighted.mean(cost_fwd_Below, w = N_Above, na.rm = T),
            `Current Event, Old` = weighted.mean(cost_current_Above, w = N_Above, na.rm = T),
            `Current Event, Young` = weighted.mean(cost_current_Below, w = N_Above, na.rm = T),
  ) %>%   
  mutate_if(.predicate = is.numeric, .funs = ~format(., digits = 0, big.mark = ",", scientific = F))



stargazer(tab_rewegihted_age_total, summary = F, out = "tab_reweighted_age_total.tex", rownames = F)



# 3d figs -----------------------------------------------------------------


# setup -------------------------------------------------------------------

floor_log10 <- function(x)
  10 ^ floor(log10(x))


library(viridis)
# import data -------------------------------------------------------------

fig3_phat0_3d <- fread("figure3_phat0_3d.csv") %>% mutate(phat_type = "Initial")
fig3_phat_dynamic_3d <- fread("figure3_phat_dynamic_3d.csv") %>% mutate(phat_type = "Current")

dt <- bind_rows(fig3_phat0_3d, fig3_phat_dynamic_3d)

# make plots --------------------------------------------------------------



render_plot <- function(p) {
  p + scale_x_continuous(breaks = seq(0, 1, .2)) +
    scale_y_continuous(breaks = seq(0, 12, 1)) +
    labs(x = paste0("Prognosis (One-Year Mortality Risk)"), y = "Month from Initial Diagnosis") +
    facet_grid(phat_type ~ group) +
    theme(panel.grid = element_blank(), legend.position = "bottom")
}

p1 <- dt %>%
  mutate(adj_cost = sum_cost / sum_days * 30) %>%
  ggplot() + geom_rect(aes(
    xmin = bins,
    xmax = bins + .1,
    ymin = time,
    ymax = time + 1,
    fill = adj_cost
  )) + scale_fill_continuous(type = "viridis")
render_plot(p1) +
  labs(fill = "Average Monthly Spending") + theme(legend.key.width = unit(1.3, "cm")) +
  ggsave(
    paste0("spending_on_phat_and_month_by_group.pdf"),
    w = 10,
    h = 7.5
  )


p2 <- dt %>%
  mutate(logN = floor_log10(N)) %>%
  ggplot() + geom_rect(aes(
    xmin = bins,
    xmax = bins + .1,
    ymin = time,
    ymax = time + 1,
    fill = as.factor(logN)
  )) +
  scale_fill_brewer(
    breaks = 10 ^ (1:4),
    labels = paste0(10 ^ (1:4), "-", c(10 ^ (2:4), "")),
    palette = "Blues"
  )
render_plot(p2)  +  labs(fill = "Number of Cases") +
  ggsave(
    paste0(
      "ncases_by_phat_and_month_by_group.pdf"
    ),
    w = 10,
    h = 7.5
  )

